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In  this  paper,  we  present  a  new  numerical  method  for  the  solution  of  linear  two- point  boundary 
value  problems  of  ordinary  differential  equations.  After  reducing  the  differential  equation  to  a 
second  kind  integral  equation,  we  discretize  the  latter  via  a  high  order  Nystrom  scheme.  A  somewhat 
involved  analytical  apparatus  is  then  constructed  which  allows  for  the  solution  of  the  discrete  system 
using  0( N  ■  p2)  operations,  where  N  is  the  number  of  nodes  on  the  interval  and  p  is  the  desired 
order  of  convergence.  Thus,  the  advantages  of  the  integral  equation  formulation  (small  condition 
number,  insensitivity  to  boundary  layers,  insensitivity  to  end-point  singularities,  etc.)  are  retained, 
while  achieving  a  computational  efficiency  previously  available  only  to  finite  difference  or  finite 
element  methods. 
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I.  Introduction 


Second  kind  integral  equations  have  been  a  popular  analytical  tool  in  the  study  of  ordinary 
differential  equations  for  nearly  a  century.  When  boundary  value  problems  are  being  considered, 
the  integral  equations  which  arise  are  of  the  Fredholm  type.  From  an  abstract  viewpoint,  the 
advantage  of  this  formulation  is  that  many  properties  of  the  solution  are  readily  apparent.  From 
a  computational  viewpoint,  the  advantage  of  this  formulation  is  that  the  linear  systems  which 
arise  from  discretization  are  generally  well-conditioned.  An  ill-behaved  differential  equation, 
such  as  a  high  order  Bessel  equation,  can  often  be  reduced  to  a  perfectly  tractable  integral 
equation  by  means  of  an  appropriate  choice  of  the  “background”  Green’s  function  (see  Example 
2  in  Section  5  below).  Standard  finite  difference  and  finite  element  methods,  on  the  other  hand, 
which  discretize  the  original  differential  equation,  encounter  serious  numerical  difficulties  when 
the  solution  possesses  derivatives  of  large  magnitude  (boundary  layers).  A  second  advantage  is 
that  there  exist  extremely  stable,  high  order  numerical  methods  for  the  solution  of  second  kind 
Fredholm  equations,  while  the  order  of  convergence  of  most  practical  schemes  for  the  solution 
of  ordinary  differential  equations  tends  to  be  limited,  even  if  Richardson  extrapolation  and 
deferred  correction  approaches  are  considered. 

Despite  all  these  advantages,  integral  equations  are  virtually  never  used  as  a  numerical  tool 
for  the  solution  of  two- point  boundary  value  problems,  since  their  discretization  leads  to  dense 
systems  of  linear  algebraic  equations,  and  the  solution  of  a  dense  linear  system  of  dimension  AT 
requires  order  0(N 3)  arithmetic  operations.  Finite  difference  and  finite  element  schemes  lead 
to  banded  systems  of  linear  algebraic  equations,  and  the  solution  of  the  latter  requires  order 
O(N)  arithmetic  operations,  where  N  is  the  dimension  of  the  problem.  This  makes  the  use  of 
integral  equations  extremely  unattractive  as  a  numerical  tool,  despite  their  superior  analytical 
properties.  A  similar  difficulty  is  encountered  when  spectral  methods  are  applied  to  boundary 
value  problems.  They  yield  high  order  accuracy,  but  result  in  dense  systems  of  linear  algebraic 
equations. 

In  [9],  it  is  observed  that  while  the  integral  operators  of  one-dimensional  potential  theorv 
are  dense,  they  can  be  applied  to  arbitrary  functions  in  a  “fast”  manner  (for  a  cost  proportional 
to  the  number  of  nodes  N).  This  observation  is  then  used  to  construct  iterative  schemes  for 
the  numerical  solution  of  second  kind  Fredholm  equations  resulting  from  two-point  boundary 
value  problems,  with  an  asymptotic  CPU  time  estimate  proportional  to  N .  Unfortunately,  the 
number  of  iterations  required  by  the  resulting  procedure  depends  on  the  problem  being  solved 
(the  usual  drawback  of  iterative  schemes),  leading  in  many  cases  to  excessive  computation 
times.  A  different  approach,  using  Chebyshev  polynomials  is  described  in  [8],  A  method  is 
constructed  which  solves  the  corresponding  integral  equation  directly.  It  requires  an  amount 
of  work  of  the  order  0(N  \ogN),  but  applies  only  when  the  differential  operator  has  constant 
coeffficients. 

In  this  paper,  we  extend  the  results  of  [9]  and  [8]  by  showing  that  not  only  the  integral 
operators  of  one-dimensional  potential  theory  but  also  their  inverses  can  be  applied  numerically 
to  arbitrary  functions  for  a  cost  proportional  to  N .  This  observation  is  used,  in  conjunction 
a  pth  order  Nystrom  scheme,  to  construct  a  fast  algorithm  for  the  solution  of  the  original 
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differential  equation.  While  the  asymptotic  CPU  time  estimate  is  proportional  to  p2.Y,  the 
algorithm  retains  the  flexibility  and  stability  expected  from  second  kind  Fredholm  equations. 

The  plan  of  this  paper  is  as  follows:  in  Section  2  we  summarize  the  relevant  properties  of 
Green's  functions  for  second  order  differential  equations,  in  Section  3  we  develop  the  analytical 
apparatus  to  be  used,  and  in  Section  4  we  describe  the  numerical  scheme  itself.  The  performance 
of  the  method  is  illustrated  in  Section  5  with  practical  examples.  Our  conclusions  and  several 
generalizations  are  discussed  in  Section  6.  Finally,  in  Appendix  A,  we  describe  the  relevant 
portions  of  the  theory  of  Chebyshev  approximation  and  quadrature. 

The  algorithm  of  this  paper  is  based  on  a  set  of  simple  observations,  but  involves  a  large 
amount  of  notation  and  a  modest  amount  of  algebraic  manipulation.  For  the  sake  of  clarity, 
we  will  attempt  to  use  two  levels  of  description  throughout,  one  cursory  and  qualitative  and 
the  other  detailed  and  rigorous. 

II.  Mathematical  and  Numerical  Preliminaries 


In  this  section  we  summarize  several  classical  results. 

2.1.  Green’s  functions  for  second  order  ordinary  differential  equations. 

We  consider  the  problem  of  determining  a  function  <p  in  C2[u,c]  which  satisfies  a  second 
order  differential  equation 

+  p(x)  •  <j>\x)  +  q(x)  •  <j>(x)  =  f(x)  (1) 

on  the  interval  [a,c]  C  R-,  where  p,q  :  (a,c)  — ►  R  are  continuous  functions,  subject  to  boundary 
conditions  of  the  form 


Cn  -<t>{a)  +  C12  •  -  £i>  (2) 

C21  •<*>(<:)  +  C22  •  <t>\c)  =  <2  •  (3) 

As  is  well-known,  the  original  equation  (1)  with  inhomogeneo-  ■“  G^undary  conditions  of  the 
form  (2),  (3)  can  be  reduced  to  one  with  homogenous  boundary  ond.tions 

C11  •  4>{a)  +  Ci 2  •  4>'(a)  =  0,  (4) 

C21  ■  +  C22  •  4>'(c)  =  0,  (5) 

by  the  addition  of  an  appropriate  linear  function  to  the  unknown  function  <p.  We  will,  therefore, 
assume  that  the  problem  to  be  solved  has  been  provided  in  the  form  (1).  (4),  (5). 

Let  us  now  consider  the  equation 

<f>"(x)  t  po(r)  •  <p'(r)  +  q0(x)  ■  4>{x)  =  0  (G) 

with  the  functions  Po,<?o  €  C5(a,c),  and  denote  by  G 0  the  Green’s  function  for  equation  (G) 
with  the  boundary  condition'.  (4)  and  (5).  The  standard  procedure  for  converting  a  two-point 
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boundary  value  problem  into  a  second  kind  integral  equation  then  consists  in  representing  the 
solution  <*>  of  the  problem  (1),  (4),  (5)  by  the  formula 


4>(x)=  f  G0(x,t)  ■  a(t)  di  (7) 

J  a 

with  o  :  [0, 1]  — ♦  R  a  new  unknown  function  to  be  determined.  Substituting  (7)  into  (1),  we 
obtain  the  desired  integral  equation 

<r(x)  +  \p(x)  ~  Po(x)}  ■  [  Gi(x,t)a(t)dt 

J  a 

+  [<?(*) -9o (*)]•/  G0(x,t)cr(t)dt  =  f(x)  (8) 

Ja 

where  the  function  Gx  :  [a,c]  x  [a,c]  —  R  is  defined  by  the  formula 

G\(x,  t)  =  ~Go(x,  t)  (9) 

for  all  (x,t)  6  [a,c]  X  [a,c]. 

Of  course,  if  po{x)  =  p(x )  and  qo{x)  =  q{x),  then  the  solution  to  equation  (8)  is  trivially 
a  =  /.  Our  working  assumption  is  that  for  some  functions  po-9o.  the  the  Green’s  function  is 
known  or  computable,  but  that  for  the  original  differential  equation  the  Green’s  function  is 
unavailable.  Fortunately,  there  is  a  well  known  mechanism  for  constructing  Green’s  functions 
from  known  independent  solutions  of  a  differential  equation.  This  construction,  described  in 
the  following  theorem,  is  the  principal  analytical  tool  of  the  paper.  For  a  proof,  see  [3]. 


Theorem  2.1  Suppose  that  ui,uT  :  [a,c]  — ►  R  are  two  linearly  independent  solutions  of  the 
equation  (6),  such  that  u\  satisfies  the  boundary  condition  (4)  and  uT  satisfies  the  condition 
(5).  Suppose  further  that  the  functions  u\,u'r,vuvT  :  [a,c]  — *  R  are  defined  by  the  formulae 


u'fit) 

<(t) 

Vl{t) 

VT{t) 


d  ,  , 

Tt*M, 


_ _ 

u'r{t)  ■  U l(t)  -  u'fit)  ■  Ur(t)' 

_ «r(Q _ 

«r(0  '«t(0  "“{(O'  MO' 


(10) 


Then  the  Green’s  function  and  its  derivative  for  equation  (1)  with  the  boundary  conditions  (4) 
and  (5)  are  given  by 


Gq(xJ)  =  ui(x)  ■  vt(t)  for  t  <  x, 

G0(x,t)  =  ur(x)-vr(t)  for  x  <  t,  (11) 


t 
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(12) 


G\[x,t)  =  u't( x)-vi(t)  for  t  <  x, 

G\(x,t)  =  u'r(x)  ■  vr(t)  for  x  <t. 

2.2.  A  Lemma  from  Linear  Algebra. 

The  following  lemma  provides  analytical  inverses  for  rank  one  perturbations  of  the  identity 
matrix.  It  is  a  particular  case  of  the  Sherman- Morrison  formula  (see,  for  example,  [6])  and  is 
easy  to  verify  directly. 

Lemma  2.1  For  any  two  vectors  {/,  V  £  L 2  such  that  (U,V)  ^  1, 

(J-U  0  VT)-*  =  I  +  --  ----  •  U  o  VT.  (13) 


III.  The  Analytical  Apparatus 

In  the  remainder  of  this  paper,  we  assume  that  the  solution  to  the  differential  equation 
(1)  is  being  sought  on  the  interval  [a,c]  and  that  b  is  some  intermediate  point  (a  <  b  <  c). 
The  fundamental  observation  on  which  the  fast  algorithm  is  based  is  that  the  solution  to  the 
integral  equation  (8)  on  the  entire  domain  [a,c]  can  easily  be  constructed  from  the  solution  of 
two  independent  integral  equations,  one  defined  on  [a,  b]  and  one  on  [6,c].  This  leads  naturally 
to  a  recursive  algorithm,  in  which  independent  solutions  on  a  large  number  of  subintervals  are 
successively  merged  until  the  full  solution  is  obtained.  A  precise  formulation  of  the  construction 
and  the  resulting  numerical  scheme  will  require  some  notation. 

3.1.  Notation. 

We  will  denote  the  subintervals  [a,  6]  and  [6,c]  of  [a,c]  by  A  and  B,  respectively.  For 
convenience,  we  write  the  integral  equation  (8)  in  the  form 

o(x)  +  p(x)-  f  Gl{x,t)cr(t)dt  +  q[x)  ■  f  G0{x,t)a{t)  dt  =  f{x) 

Ja  J  a 

where  p(i)  =  p(x)  -  po(x)  and  q(x)  =  q(x)  -  qo(x).  The  functions  Go,Ci  .’  [a,c]  X  [a,c]  — *  R 
are  the  Green's  function  and  its  derivative  defined  by  formulae  (11)  and  (12). 

We  define  the  operator  P  :  L2[a,c ]  — ►  L2[a,c]  by 

P(a)(x)  =  o(x)  +  p(x)  ■  f  G-[(x,t)  ■  cr(<)  dt  +  q(x)  ■  f  G0(x,t)  ■  <r(t)  dt  ,  (14) 

J  a  Ja 

so  that  equation  (8)  assumes  the  form 

Pa  =  f.  (15) 

We  will  require  the  four  operators 


Pa  a  :  L2[a,  b]  — >  L2[a,  6]  , 
Pab  :  P[b,c}^  L2[a,b], 
Pba  :  L2M]-*L2[6,c], 
Pbb  ■  L2[b,c]  — >  L2[6,c] 
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defined  bv 


Paa(v)(x)  =  v{x)  +  P{?)-f  G1(x,i)-o(t)dt  +  q[x)-(  G0(x ,t)  ■  cr{t)  dt  , 

J  a  J  a 

Pab(0)(x)  =  p{x)  ■  Jb  G^(x,t)  ■  o(t)dt  +  q(x)  ■  G0(x,t)  ■  cr(t)dt  , 

Pba(v){x)  =  p{*)-f  Gj(x,0  ■  <r{t)dt  +  q{x)  ■  f  G0(x,t)  ■  a{t)dt  , 

Pbb(c){*)  =  <r(x)  +  p{x)  ■  f  Gi(x,t)  ■  a{t)dt  +  q{x)  ■  f  G0(x,  t)  ■  a(t)  dt. 

Jl  Jb 

We  will  also  require  the  functions  pi  and  rpT  defined  by 

M?)  =  p(2-)  •  u',(x)  +  $(x)  •  Ul(x) 

il>r(x)  =  p(x)  ■  u'T(x)  +  q(x)  ■  ur(x)  . 

Given  a  function  /  €  L2[a,c],  we  will  follow  the  convention  of  denoting  its  restriction  to  A  and 
B  by  f\A  and  /|g  respectively.  Assuming  that  the  operators  P,Paa,Pbb  are  non-singular,  we 
then  define  the  mappings 


(16) 

(IT) 

(IS) 

(19) 

(20) 
(21) 


via  the  formulae 


XI -\r 

:  [a,c]  —  R, 

</>lA,d>rA 

:  A  —  R  , 

.  B  -*  R 

XI  = 

p-'m, 

(22) 

Xr  “ 

P-'Wr), 

(23) 

K  = 

Paa^'I\a)^ 

<t>rA  = 

Paa  ( ^’r\ a  )’ 

<filB  = 

Fgg(l/’l|B), 

(24) 

<t>rB  = 

Fgg(t/-V|B)-  • 

Finally,  we  will  define  three  2x2  matrices  aA,oB  and  a  by  the  formulae 

Qn  =  =  {viIA,<t>rA)  , 

Q21  ={Vr,A,<PlA)  ,  «22  =  (^  >A,<t>rA), 


(25) 


Ol  2  =  {VliB,<}>TB)  , 
Q22  =  (Vr\ Bi<t>TB)  i 


ofi  =  (viIB,<t>iB)  , 
of\  =  {vT\B,4>lB)  , 
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(26) 


an  =  (l’/,X/)  ,  ai2  =  (vt,\r), 

«21  =  (»->.  XI )  i  Q22  =  (l’r,Xr).  (27) 

and  the  coefficients  St.ST  by 

61  =  (vi,a), 

6r  =  (  vr,a),  (28) 

where  <r  is  the  solution  to  equation  (15). 

3.2.  Analysis  of  the  operators  Pab- Pba- 

In  this  subsection,  we  observe  that  each  of  the  operators  Pab  an<3  Pba  is  °f  rank  one  and 
give  simple  expressions  for  these  operators  as  outer  products. 


Lemma  3.1  In  the  notation  of  the  preceding  subsection, 

Pab  =  Vv|yl  o  vJ[B,  (29) 

Pba  =  i’t\B  0  vIlA-  (3°) 

Proof.  Combining  expressions  (11),  (12)  and  (17),  and  observing  that  x  <  t  for  any 
x  £  [a,  6],  t  £  [6,  c],  we  have  for  any  o  £  L2[b ,  c], 

Pab(<t)(x)  =  p(x)  ■  Gi(x,t)  ■  o(t)  dt  +  ?(i)  ■  J  Go{x,t)  •  <r(t)  dt 

=  p(x)-  J  u'r(x)  ■  vr(t)  ■  o(t)  dt  +  g(x)  ■  uT(x)  ■  vT(t)  ■  cr(t)  dt 

=  (p(x)-u' r(x)  +  g(x)-ur(x))-(vriB,<r).  (31) 

The  result  (29)  now  follows  from  the  definition  of  ipr  in  equation  (21). 

Similarly,  combining  expressions  (11),  (12)  and  (18),  we  observe  that  for  any  cr  €  L2[b,c) 
and  x  £  [h,  c], 

Pba(°)(x)  =  P(x)  ■  J  Gi(x,t)  ■  o(t)  dt  +  <j(x)  ■  J  G0(x,t)  ■  o(t)  dt 

=  p(x)  ■  J  u't(x)  ■  Vi{t)  ■  o(t)  dt  +  q(x)  ■  J  u,{x)  ■  v,{t)  ■  o(t)  dt 

=  iP(x)  ■  «{(*)  +  q(x)  ■  ui(x))  •  (ni|X,CT).  (32) 

The  result  (30)  follows  from  the  definition  of  rfi  in  equation  (20).  D 
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3.3.  Recursive  solution  of  the  integral  equation  (15) 
We  now  consider  the  original  integral  equation  (15) 

Pa  =  f 


and  two  auxilliary  equations 

PaakVa)  =  f\A  (33) 

Pbb(^b)  =  f[B  ■  (34) 

The  main  result  of  this  subsection  is  the  following  lemma,  which  constructs  the  solution  a 
of  equation  (15)  from  the  solutions  t]a,VB  of  equations  (33)  and  (34). 


Lemma  3.2  If,  in  the  notation  of  Subsection  3.1,  all  three  operators  P,  Paa-,  Pbb  are  non¬ 
singular,  then 

a\A  =  Va~  (6?  ~  ~  (if  -  a?2-S?))-<J>TA,  (35) 

<>\B  =  V B~  ~  ^  •  (6?  -  a?,  •  6f))  ■  4>lB ,  (36) 

where  the  real  numbers  Sf,  A®,  6f,  bf,  A  are  given  by 

tf  =  0vw)’ 

6r  =  (%4,*M),  (37) 

=  (»,|B,f®), 

A  =  1  -of,  •«£.  (38) 

Proof.  Using  definitions  (14)  -  (19),  the  integral  equation 

P°  =  f 


can  be  rewritten  in  the  form 


Paa{<?\a)  +  Pab{o\B)  =  f\A,  (39) 

Pba(<*\a)  +  Pbb(p\b)  =  f\B-  (40) 

The  outer  product  expansions  (29)  and  (30)  for  Pab  and  Pba •  respectively,  can  then  be  used 
to  obtain  an  explicit  solution  to  the  coupled  equations  (39)  and  (40)  in  terms  of  the  functions 

Va<  Vb<  4>ia >  4>ta,  4>ib  and  <f>TB.  Indeed,  applying  the  operator  P^\  to  equation  (39)  and  the 

operator  PBB  to  equation  (40),  we  have 

a\ A  +  PaA  0  Pab{<?\b)  =  Paa(/\A  )>  (41 ) 
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(42) 


PBB  0  Pba(C\A)  +  <T| B  =  p bb(/|s). 

Substituting  the  outer  product  expansions  (29)  and  (30)  into  (41)  and  (42)  yields 


”IA  +  Pa  a  0  VVM  o  vJ[B  o  alB  =  va, 

PBB°  *l\y  0  0  °\A  +  a\B  =  VB'  ( 41 ' 

or 

0\A  +  <t>rA  °  vJiB  0  o\B  -  r)A,  (4'. ) 

4>ib  O  vfiA  0  crjA  +  OfB  =  i IB,  (4(; ) 

where  we  have  used  the  definitions  (24)  for  <prA  and  4nB.  Now,  multiplying  (46)  by  <pTA  o  r J 
and  subtracting  it  from  (45),  we  obtain 

(I  ~  <t>rA  0  vjxg  o  <t>iB  o  )  a{A  =  r?4  -  <pTA  o  0  VB-  ( 47 ) 

Similarly,  multiplying  (45)  by  d>;B  o  and  subtracting  it  from  (46),  we  get 

(7  -  <t>,B  o  o  <j>rA  o  t^B)  <T|B  -r)B-  <p‘B  o  f/^  o  T)A .  (48) 

Due  to  (25)  and  (26),  we  can  rewrite  these  equations  in  the  form 

V  -  ofj  •  </>r„  O  )  ^|.4  =  °  V*B  0  T?B,  (49  ) 

(7  -  af2  •  <ftB  o  v*A )  o\B  ~m~  <t>iB  0  VT[A  0  VA-  (5°) 

By  application  of  Lemma  2.1,  we  obtain 

aB 

°\A  =  (/  +  •  <t>rA  0  vflA)  (r)A  -  <j>rA  o  Ur)fi  O  T?S),  (51) 

^|S  =  (/  +  ^  •  <PiB  0  ^B)  (r?B  -  <t>lB  O  O  tm).  (52) 

The  results  (35)  and  (36)  now  follow  from  simple  algebraic  manipulation  using  equations  (25). 
(26)  and  (37).  □ 


Remark  3.1.  Suppose  that  6i  and  b2  are  a  pair  of  real  numbers  such  that  a  <  bj  <  b2  < 
c,  and  that  the  interval  [6^  62]  is  denoted  by  C.  We  will  denote  by  Pec  the  restriction  of 
the  operator  P  to  the  interval  C.  Assuming  that  Pec  is  invertible,  we  define  the  functions 
Vc,<t>ic,4>rc  ■  C  — *  R  by 


Vc 

l8 

II 

(53) 

<t>ic 

—  ^Cc(^fic)  » 

(54) 

tre 

=  PCC  ( c )  • 

(55) 
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By  applying  the  above  lemma  twice  (once  for  the  subinterval  [a, 61]  and  once  for  [0,62]).  we 
may  easily  observe  that  there  exist  two  real  numbers  Aj,  A2  such  that 

0(x)  =  Vc(x)  +  A!  -  4>ic(x)  +  A2  •  4>rc{x)  (56) 

for  all  1  £  C.  The  exact  expressions  for  the  coefficients  Ai,A2  are  complicated,  but  irrelevant 
for  the  purposes  of  this  paper.  The  existence  of  a  relation  of  the  form  (56).  however,  will  be 
critically  important  in  Section  4. 

3.4.  Further  Analytical  Results 

We  now  collect  a  number  of  identities  which  are  necessary  for  the  algorithm  to  be  presented 
in  Section  4.  First,  we  apply  Lemma  3.2  to  the  particular  cases  /  =  V’/,  /  =  to  obtain 
analytical  expressions  for  the  functions  \i  and  yr  defined  in  equations  (22)  and  (23). 

Corollary  3.1  If,  in  the  notation  of  Subsection  3.1,  all  three  operators  P,  P^a,  Pbb  are  non¬ 
singular ,  then 


Ah* 

=  4>Ia  ■ 

_(±Z 

-  Qn) 
A 

•  aB 
Q21 

■<PrA. 

(57) 

Afi  B 

=  (1-- 

°n  - 

-  ' 

A 

-*> 

■  <t>lB  ~ 

1  _ 

A 

(58) 

=  (1- 

Q22  ~ 

-  al2  ■ 

A 

nB 

—11) 

■  ^  = 

1  -  Q22 
A 

■*rA, 

(59) 

_(!: 

-  0^2) 
A 

•<*12 

'4>lg  * 

(60) 

where  the  coefficients  a *  and  a®  are  given  by  equations  (25)  and  (26). 

We  will  also  require  analytical  expressions  for  the  inner  products  61  and  6r  defined  in  (28) 
in  terms  of  the  restricted  inner  products  ,  bf  and  6?  defined  in  (37). 


Corollary  3.2  If,  in  the  notation  of  Subsection  3.1,  all  three  operators  P,  Paa.  Pbb  are  non¬ 
singular,  then 


<5,  =  (t;/,£T)=  (vi{a,<T[A)  +  (vi[b,C7IB) 

1  ~  Qfl  cA  I  S.B  ,  (0fl-l)’a12  iB 
_  +6l  + - - - 6r  . 


(61) 


6r  =  (vt,<T)  =  (tV,*,^)  +  (VTW,0\B) 

1  -  C»22  cB  ,  t.A  ,  (q22  ~  1)  '  Qfl  cA 

-  '6r  +  K  +  - £ - Si 


(62) 


9 


Proof.  Multiplying  equation  (35)  by  rj  and  tv|A -  and  equation  (36)  by  r/|g  and  tr|B,  w< 
obtain 


1 


(r'|e-^|B)  =  if  ~ 


y  B  nB  ,  r.A 

bl  iA  ,  Q11  °12  tB 

T  '  ci  +  r~  '  °r  ’ 


A 

/v  >4  4  /v  B 

yB 


(63) 

(64) 

(65) 

(66) 


A  r  A 

Adding  up  the  first  pair  of  equations  (63)  and  (64),  we  obtain  the  result  (61).  Adding  up 
the  second  pair  of  equations  (65)  and  (66),  we  obtain  (62).  D 

A  special  case  of  Corollary  3.2  is  obtained  when  /  =  tl'1  or  /  =  Vv-  The  proof  follows  easily 
from  the  definitions  of  \-/  and  \r  in  (22)  and  (23). 


Corollary  3.3 

singular,  then 


If,  in  the  noiation  of  Subsection  3.1,  all  three  operators 


,  .  x  (1  -  ofi)  -  (aft  -  Qn-Q2\'l  ,  „ b 

On  =  (tV,  \l)  =  - ^ -  +  Qn- 

,  1  Q?i  ■  ( 1  ”  °22)  •  ( 1  “  Qn  )  ,  A 

021  =  (tv,  XI )  =  - ^ -  +  Q2l- 

/  \  0^2  •  (1  -  af2)  •  (1  -  afi)  B 

«12  =  (t’l,  Xr)  =  - ^ -  +  Ql2, 

(1  ~  022)  •  (q®2  -  0^2  •  of]  ) 

022  =  (VI,  Xr)  =  - T -  +  °22- 


P.Paa.Pbb  ore  non- 


(67) 

(6S) 

(69) 

(70) 


Finally,  combining  Lemma  3.2  with  the  expressions  (57)-(60),  we  have 

Corollary  3.4  Suppose  that  in  the  notation  of  Subsection  3.1,  all  three  operators  P.  Paa  •  Pbb 
are  non-singular.  Suppose  further  that  the  function  F  :  [a,c]  — *  R  is  defintd  by  the  formula 

F{x )  =  Ai  •  x,(t)  +  A2  -Xr(x)  +  A3  •  <t(j).  (71) 


Then  on  the  interval  [a,  b], 

F(i)  =  pi  ■  <t>iA(x)  +  P2  ■  <t>rA{x)  +  M3  ■  nA{x).  (72) 

with  the  coefficients  pi ,  P2 .  pz  defined  by  th>  formulae 


Ml 

= 

Ai, 

M2 

= 

_(Jj 

-  <»n) 
A 

•o& 

■  Aj 

1  -  a22 
A 

+ 

/  °21 

1  A 

-(bf- 

-  af2 

■6*. 

)-«f)-A3 

M3 

= 

a3. 

(73) 
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(74) 


(75) 


IV.  Description  of  the  Algorithm 

We  turn  now  to  the  construction  of  the  fast  algorithm  for  the  solution  of  the  integral 
equation  (15) 

Po  =  /  , 

based  on  the  apparatus  developed  in  Section  III.  The  main  too]  at  our  disposal  is  the  ability 
to  merge  the  solutions  of  restricted  versions  of  the  integral  equation  in  adjacent  subintervals 
(Lemma  3.2).  As  this  suggests  a  recursive  procedure,  we  begin  by  subdividing  the  whole 
interval  [a,c],  on  which  the  solution  to  (15)  is  sought,  into  a  large  number  of  subintervals.  For 
the  sake  of  simplicity,  we  assume  that  m  is  a  positive  integer  and  that  M  =  2m  is  the  number 
of  subintervals  created.  The  boundary  points  of  the  subintervals  are  then  given  bv  a  strictly 
increasing  sequence  of  numbers 

b\.b2,  •  •  •  > b\i,  6a/ +i  (76) 

with  bx  =  a  and  6m+i  =  c.  We  define 

B?  =  [^n&.+i]  for  i  =  1, . . .,  M  (77) 

and  create  a  hierarchy  of  intervals  by  recursively  merging  adjacent  pairs.  That  is,  for 
k  =  m  -  1, . . .,  1.0,  we  define 

Bf  =  t\  (J  S2\+1  for  «  =  1, . . . ,  2*  .  (78) 

We  will  refer  to  each  '.  •",d  A;  as  a  level.  We  will  also  refer  to  the  two  fine  intervals  52,+ \  and 
Bi*'  as  children  ar  x.  the  larger  interval  Bf  as  a  parent. 

It  is  obvious  that 

and  that  for  each  level  k, 

2k 

!a,c]=U5f  ■  (80) 

«=i 


Similarly,  on  the  interval  [6.c], 

F(x)  =  vx  ■  <t>iB(x)  +  v2  ■  4> rB(x)  +  ■  r?s(i), 

with  the  coefficients  ux,  v2,  V3  defined  by  the  formulae 


,.  _  0  °ii)  \  ,  (1  _  <*22)  ‘  Qi42  x 

~ - ^ - +  - £ - A2 


a 


+  A3, 


v2  =  A2, 
^3  =  A3. 
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4.1.  Notation. 

Generalizing  the  notation  of  Section  III,  we  will  denote  by  P,,k  the  restriction  to  the  interval 
Bk  of  the  integral  operator  P,  so  that 


PiJ c(^K^) 


P 14,  2m_k 

=  <r(x)  +  p(x)-  Gi(x,t)  ■  <j(t)  dt 

+  t  2m-fc 

Go(x,t)  ■  cr(t)  dt 

14(.-1)  2m~k 


(81) 


for  any  <7  6  L7(Bk)-  For  each  B,  we  will  define  the  functions  77,^,  d>it  k,  d>Tt  k  :  Bk  —  R  as  the 
solutions  of  the  equations 

Pa(va))  =  V-  (82) 

PiA<t>i<.k)  =  (83) 

Pi,k(<i>T,  k)  —  .  (S4) 

provided  that  the  operator  P,  k  is  non-singular. 

Remark  4.1.  Suppose  now  that  the  operator  P,k  is  non-singular  on  the  interval  Bk .  Then, 
by  equation  (56),  there  exist  two  numbers  Aj’fc,  X*  £  R  such  that 


A*)  =  v,A*)  +  •  <&,,*(*)  +  Ay*  ‘  <Pr,Ax) 

for  all  x  €  Bk . 

For  each  k  =  0,1 , . . . ,  m,  and  i  =  1, 2, ...  ,2*,  we  define  a  2  X  2  matrix  a'’k,  by 

Qif  =  (vliBk’<h,.k)' 

°2\  = 

‘  * 

<*n  = 

a22  =  (Vr.Bki<Pr,  k), 

I 

and  the  vector  6l  k  =  (6j'k,  <>*’*)  by 

=  (vi  k.V,,k)> 

1 

<fr  =  (Vrigk  ’  V'A- 


(85) 


(86) 


(87) 


4.1.  Discretization  of  the  Restricted  Integral  Equations. 

Choosing  an  integer  p  >  1,  we  construct  the  p  scaled  Chebyshev  nodes 


r,  = 


cos 


(2;-l),]  +  (^)  .oli2 . p 


2? 


(88) 
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on  each  of  the  intervals  B”' ,  i  =  1,2,...,  M.  We  then  discretize  the  three  integral  equations 
(82).  (83)  and  (84)  via  a  Nvstrom  algorithm  based  on  p-point  Chebyshev  quadrature  (see 
Appendix  A).  The  resulting  approximations  to  the  functions  p,,*,  4>i,  k*  <£r,  *  at  the  nodes  rJt 
will  be  denoted  bv 


fh.k 

respectively. 

Remark  4.2.  It  is  well-known  that  the  order  of  convergence  of  the  approximations 
Pi.fc.  <J>|,  k,4>r,  k  to  the  functions  T)t  k,<t>i,  k,<t> r,  k  is  P-  Since  all  subsequent  steps  in  the  construction 
of  an  approximate  solution  a  to  the  integral  equation  (15)  are  analytic,  the  convergence  rate 
of  the  full  algorithm  depends  entirely  on  the  parameter  p.  For  example,  by  using  16  scaled 
Chebyshev  points  on  each  subinterval  at  the  finest  level,  one  obtains  a  sixteenth  order  method. 

4.2.  Informal  description  of  the  algorithm. 

We  begin  by  directly  solving  the  three  integral  equations  (82),  (83)  and  (84)  on  each  subin¬ 
terval  5,m  at  the  finest  level,  as  discussed  in  the  preceding  subsection.  Equation  (85)  then 
shows  that  a  restricted  to  B™  can  be  expressed  as  a  linear  combination  of  the  three  solutions 
Oi,  m .  <t>r,  m-  Thus,  it  remains  only  to  determine  the  two  coefficients 

\t,m  \«,m 
Aj  ,  A2 

for  each  of  the  M  subintervals  B™ .  Fortunately,  this  can  be  done  recursively.  To  see  this, 
suppose  that,  at  some  coarse  level  k  <  m  -  1,  we  are  given  the  coefficients  for  the 

subinterval  Bk .  Then  Corollary  3.4  provides  formulas  for  the  calculation  of  the  corresponding 
coefficients 

A2i-i.fc+i,A»-i.*+i  and  A*'*+\A*’fc+l 
for  the  two  child  intervals  Bk ^  and  ££,+1,  respectively.  For  initialization,  observe  that 


A}’0  =  =  0  (89) 

(i.e.  the  solution  of  equation  (82)  on  the  whole  interval  [a,c]  is  simply  a). 

In  order  to  use  the  formulas  (73  and  75)  of  Corollary  3.4,  however,  we  need  the  matrices 
a2t-i,fc+i^Q2i,fc+i  and  vectors  ^2f-i,fc+i,^2i,fc+i.  These  quantities  are  also  computed  recur¬ 
sively  but  in  the  opposite  direction,  namely,  from  the  finest  level  to  the  coarsest.  They  are 
certainly  available  at  level  m  directly  from  the  definition  (86).  For  the  interval  B f  at  any 
coarser  level  k  <  m  —  1,  Corollaries  3.2  and  3.3  describe  how  c*1'*  and  Sl,k  are  obtained  from 
the  matrices  a  and  vectors  <5  of  the  two  child  intervals. 

To  summarize,  the  algorithm  consists  of  three  parts.  First,  a  sufficiently  fine  subdivision 
6i,  62, . . . ,  6m+i  of  the  interval  [a,  c]  is  chosen  so  that,  on  each  of  the  intervals  £,,m,  the  functions 
Vt,m  1  m ,  and  d>r,,m  can  be  accurately  represented  by  a  low  order  Chebys,  -  v  expansion.  On 


each  of  the  intervals  B,_ m,  the  equations  (82)  -  (84)  are  solved  (approximately)  by  direct 
inversion  of  the  linear  system  arising  from  a  Nystrom  discretization.  Second,  the  matrices 
a'  k  and  vectors  S'k  are  computed  in  an  upward  sweep,  beginning  at  the  finest  level  m.  Finally, 
the  coefficients  \\’k  and  are  computed  in  a  downward  sweep,  beginning  at  the  coarsest 
level.  The  desired  function  a  is  then  recovered  on  each  subinterval  from  equation  (85). 

The  following  is  a  more  detailed  description  of  the  numerical  procedure. 

Algorithm 

Comment  [Define  computational  grid  ] 

Create  M  =  2m  subintervals  on  [a,c]  by  choosing  a  sequence  of  boundary  points  b^b?,. .  .  ,1>ju  i.'t  +  i 
with  bi  =  a  and  fcw  +  i  =  c  Choose  the  number  p  of  Chebyshev  nodes  on  each  interval  B™  —  [6,.  61  +  i] 
for  i  =  1, . . . ,  M .  Determine  the  locations  of  the  scaled  Chebyshev  nodes  rl ,  r2, .  .  . ,  r,p  on  each  interval 

sr . 


Step  1. 

Comment  [Construct  the  approximate  solutions  m ,  4>r<  m  of  equations  (82)  *  (84)  on  each  in¬ 

terval  B,m.] 

do  i  =  1,2 . M 

(1)  Construct  the  three  p  x  p  linear  systems  on  Bk  obtained  through  a  Nystrom 
discretization  of  the  corresponding  integral  equation  (see  Appendix  A,  Section  3). 

(2)  Solve  the  three  p  x  p  linear  systems  on  Bk  by  Gaussian  elimination, 
obtaining  the  values  r)l  m , 

enddo 


Step  2. 

Comment  [Construct  the  coefficients  a]'”,  a\'™ ,  a'^™,  6[m,  6J.,m  on  each  interval  at  the  finest 

level] 

do  i  =  1,2....,  A/ 

Evaluate  the  coefficients  a]  ™,  a^'™, ^r'm  by  applying  the  k— point 
Chebyshev  quadrature  formula  (Appendix  A,  Section  2)  to  the  inner  product  integrals  (86),  (87). 
end  do 


Step  3  (Upward  Sweep). 

Comment  [Construct  the  matrices  a'k  and  the  vectors  S' 1  for  all  intervals  at  all  coarser  levels  k  = 

m  -  1,  m  —  2, . . . ,  0.] 

do  k=  m-1,  0,  -1 
do  i=l,  2^ 

Use  formulae  (67)  -  (70)  and  (61)  -  (62)  to  compute  the  matrix  a'k 
and  the  vector  6,  k  from  the  corresponding  data  in  the  two  child 
intervals  (Qr2’~J'i+1,a2*>i+,,62’~1’l+1,62’-*+l). 

end  do 
end  do 


14 


Step  4  (Downward  Sweep). 


Comment  [Construct  the  coefficients  A,1'm,  Aj™  for  all  intervals  at  the  finest  level.) 

Set  A®'1  =  A®1  =0  (see  equation  (89)  above). 

do  k=0,m-l 
do  i=l,  2^ 

Use  Corollary  3.4  to  compute  the  coefficients  Aj+I’2  *-1,  A*+1'2  ,_I,  Aj,,ir+1,  A2*'*+1 
for  the  child  intervals  and  from  the  coefficients  A\,fc ,  Aj*  of  the  parent  interval  B*. 

end  do 
end  do 


Step  5. 

Comment  [Compute  the  solution  a  of  equation  (15)  at  the  nodes  r,1,  r-2, . . .,  rp  for  each  interval  By 
at  the  finest  level.) 

do  i=l,  M 
do  j=l,p 

Determine  the  values  of  the  solution  a  of  equation  (15)  at  the  node  rf  via  formula  (85). 
end  do 
end  do 


Step  6. 

Comment  [Compute  the  solution  <t>  of  equation  (1)  and  its  derivative  <j>'  from  the  values  of  <r.) 

Evaluate  the  integral  (7),  and  its  derivative  by  using  composite  Chebyshev  quadrature 
(see  Remark  4.4  below). 

Remark  4.3.  Inspection  of  the  above  algorithm  shows  that  the  amount  of  work  required  is 
of  the  order  0(M  p3).  Three  p  x  p  linear  systems  have  to  be  solved  for  each  of  the  M  intervals 
B?1  in  Step  1,  while  Steps  2  -  5  require  no  more  than  0(M  •  p 2)  operations.  Since  N  =  M  •  p 
is  the  total  number  of  nodes  in  the  discretization  of  the  interval  [ a,c ],  we  can  write  the  CPU 
time  estimate  in  the  form  0(N  •  p2).  The  cost  of  evaluating  the  solution  <f>  of  the  differential 
equation  (1)  from  the  integral  representation  (7)  is  (^(^logp)  (see  Remark  4.4  below). 

Remark  4.4.  The  final  step  in  the  algorithm  involves  the  evaluation  of  integrals  of  the 
form  (7)  at  each  of  the  Chebyshev  nodes  t\  on  each  subinterval  By ,  namely 


4>{r3x)  = 

[  Go{T-,t)  ■  <7(0  dt 

J  a 

(90) 

<W)  = 

f  Gi(rf  ,t)  •  <r(t)  dt  . 

Ja 

(91) 
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If  these  integrals  were  calculated  independently  for  each  t-  ,  the  amount  of  work  required 
would  be  of  the  order  0(A'2),  and  would  dominate  the  construction  of  the  function  a.  In  fact, 
this  is  unnecessary,  for  we  may  write 


4.  r1, 

vi(t)  ■  a(t)di  +  /  vi(t)  ■  a{i)  dt 
J  b, 

+«r(r/)  • 

[+/,  vT(t)  ■  <r(i)  dt  +  vr(t ) 

=  «{(*?)•[/ 
Ja 

b<  r! 

v\(t)  ■  cr(t)  dt  +  J  vi(t)  ■  a(t)  dt 

+  f  vT(t)  ■  cr(t)  dt  +  /  vT{t) 

Jrj  Jb,+ 1 

(92) 


(93) 

(9-1) 


where  we  have  used  the  representations  (11)  and  (12)  and  the  fact  that  rt7  lies  in  the  interval 
=  [6,.6,+1].  Step  6  can  then  be  written  in  detail  as  follows: 


Step  6  (a). 

Comment  [Precompute  the  integrals  of  vra  and  vr  a  on  each  subinterval  S,m  by  Chebyshev  quadrature. 
These  integrals  will  be  denoted  U  and  Ir ,  respectively  ] 

do  i=l,  M 

h{B?)  =  /‘'+1  v,(t)a(t)dt. 

end  do 


Step  6  (b). 

Comment  [March  across  interval  from  a  to  c,  computing  <f>  and  <j>'  at  each  node  in  discretization.  The 
varaibles  J\  and  3r  will  be  used  to  accumulate  the  integrals  J*'  vi(t)  ■  <?(t)dt  and  j  tv(<)  -  cr(t)dt, 
respectively. 

Set  Jr  =E"2Ir(Br). 

Set  3  =  0. 

do  i=l,  M 
do  j=l,p 

For  each  rf ,  compute 

<HH )  =  «f(^)  •  \ji  +  //;  V,(t)  ■  ff(t)  dt]  +  ur(rf )  ■  |/t4;+’  Vr{t)  ■  a {t)dt  +  Jr  j 

*(H)  =  [?•  +  Mt)  *(<)*]  +  W)  ■  [/t6;+’  MO  •  *(0 dt  +  Jr] 

end  do 
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Ji  =  Ji  +  h(B™ ) 
Jr  =  Jr~  Ir(B?+i) 
end  do 


Thus,  the  amount  of  work  required  in  Step  6(a)  is  O(N).  The  integrals  required  on  each 
subinterval  in  Step  6  (bl  can  be  computed  by  spectral  integration  (see  Appendix  A,  Section  2) 
using  O(plogp)  work.  The  total  cost  is  therefore  of  the  order  O(Mplogp)  or  0(Arlogp). 


V.  Numerical  Results 

FORTRAN  programs  have  been  written  implementing  the  algorithm  described  in  the  pre¬ 
ceding  section,  for  both  real  and  complex  valued  functions.  In  this  section,  we  discuss  several 
details  of  our  implementation,  and  demonstrate  the  peformance  of  the  scheme  with  four  nu¬ 
merical  examples. 

The  following  technical  details  of  our  implementation  appear  to  be  worth  mentioning. 

1.  Originally,  the  algorithm  was  implemented  for  a  fairly  wide  choice  of  background  equations 
(6),  and  corresponding  Green’s  functions  (11).  Our  numerical  experiments  showed  that  the 
advantages  of  one  background  Green’s  function  over  another  tend  to  be  minor,  unless  the 
original  equation  (1)  and  the  equation  (6)  from  which  the  background  Green’s  function  is 
constructed  can  be  chosen  to  be  extremely  close.  Therefore,  all  subsequent  implementations 
and  all  numerical  experiments  were  performed  with  equation  (6)  of  the  form 

4>"(x)  =  0.  (95) 

2.  The  algorithm  described  in  the  preceding  section  requires  that  the  number  M  of  elementary 
subintervals  on  the  interval  [a,  c]  be  a  power  of  2.  Cleary,  this  is  not  an  essential  limitation 
and  it  can  be  removed  by  simple  bookkeeping  changes.  In  the  version  of  the  algorithm  used 
for  numerical  experiments,  these  changes  were  made. 

3.  The  algorithm  depends  for  its  stability  on  the  equations  (82)  -  (84)  having  unique  solutions 
for  all  subintervals  Bf  ( k  -  0,1,..., M,  t  =  1,...,2*.)  It  is  easy  to  construct  examples  for 
which  this  condition  is  violated,  even  though  equation  (15)  has  a  unique  solution.  In  such  cases, 
a  different  subdivision  of  the  interval  [a,c]  can  be  attempted,  such  that  none  of  the  subintervals 

of  the  new  subdivision  coincides  with  an  interval  of  the  original  one.  This  procedure  can  be 
viewed  as  a  form  of  pivoting,  and  it  is  easy  to  show  that  it  is  always  possible  to  make  it  work. 
It  has  not  been  implemented  at  this  point,  and  we  have  not  so  far  encountered  a  need  for  it. 

4.  We  have,  however,  implemented  a  crude  scheme  for  detecting  high  condition  numbers  in 
the  algorithm.  These  can  occure  in  two  places:  in  the  solution  of  the  linear  systems  on  each 
of  the  finest  level  subintervals  (Step  1),  and  while  merging  the  solutions  on  two  consecutive 
subintervals  via  formulae  (61)-(62)  and  (67)-  (70)  (Step  3).  In  the  first  case,  the  condition 
number  of  the  system  being  solved  is  estimated  in  the  process  of  solution  (we  use  a  standard 
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UNPACK  routine),  and  the  largest  of  these  is  returned  to  the  user.  In  the  second  case,  the 
immediate  reason  for  the  ill-conditioning  is  the  appearance  of  small  values  of  the  coefficient  A 
in  formulae  (61)-(62)  and  (67)-(70).  The  smallest  of  these  is  also  returned  to  the  user.  When 
an  extremely  large  condition  number  is  detected  by  the  UNPACK  routine,  or  an  extremely 
small  A  is  generated  in  the  merging  process,  the  resulting  solution  of  the  original  ODE  should 
be  viewed  as  suspect.  It  is  easy  to  show  that  when  the  differential  operator  is  positive  definite, 
this  cannot  happen.  A  more  complete  treatment  of  this  subject  requires  further  study. 

5.  In  the  upward  sweep  (Step  3),  we  evaluate  the  matrices  u,fc  for  all  intervals  Btik  and  use 
these  matrices  to  evaluate  the  vectors  6,  the  coefficients  A,  and,  finally,  the  solution  a  of  the 
integral  equation  ( 15).  But  the  matrices  a,,k  do  not  depend  on  the  right-hand  side  /  of  equation 
(15),  and  it  is  easy  to  see  that  their  evaluation  accounts  for  more  than  three  quarters  of  the 
work.  Therefore,  whenever  the  equation  (15)  has  to  be  solved  with  multiple  right-hand  sides, 
we  precompute  the  matrices  a*’*  and  store  them,  saving  about  75%  of  the  cost  of  the  evaluation 
of  subsequent  solutions. 

The  algorithm  of  this  paper  has  been  applied  to  a  variety  of  problems.  Four  experiments 
are  described  below,  and  their  results  are  summarized  in  Tables  1-10.  In  each  of  these  tables, 
the  first  column  contains  the  total  number  JV  of  nodes  in  the  discretization  of  the  interval 
[a,c].  The  second  column  contains  the  relative  L 2  error  of  the  numerical  solution  as  compared 
with  the  analytically  obtained  one,  and  the  third  column  contains  the  maximum  absolute  error 
obtained  at  any  node  in  the  discretization.  Columns  four  and  five  contain  the  same  information 
for  the  derivative  of  the  solution  (i.e.  its  relative  L 2  and  /,«>  errors  respectively).  Finally,  the 
last  column  contains  the  CPU  time  required  to  solve  the  problem.  In  all  cases,  the  times  given 
are  for  a  SUN  3/60  computer  using  the  68881  floating  point  coprocessor. 

Example  1.  This  example  is  taken  from  [10],  where  it  is  described  as  a  reasonably  difficult 
one  due  to  the  presence  of  rapidly  growing  solutions  of  the  corresponding  homogenous  equation. 
The  equation  to  be  solved  is 

<t>"  +  400d>  =  -400cos2(7rx)  -  2n2  cos(2x  x)  (96) 

with  the  boundary  conditions 

=  4>(1)  =  0.  (97) 

The  algorithm  has  been  applied  to  this  problem  with  p  =  8, 16  and  24,  and  the  results  of  this 
experiment  are  presented  in  Tables  1-3. 
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n 

wm  hi  ii 

£°°(0) 

EV) 

t  (sec.) 

16 

0.954  Xl0~5 

0.659  Xl0-5 

0.959  XlO"5 

0.138  X10“3 

0.200  xl0° 

32 

0.457  xlO'8 

0.301  XlO-8 

0.545  X10~8 

0.602  XlO"7 

0.340  xl0° 

64 

0.401  XlO-12 

0.388  XlO-12 

0.581  xl0~12 

0.778  xl0~n 

0.700  x  10° 

128 

0.658  X10-15 

0.139  x  10-14 

0.106  XlO-14 

0.319  xl0~13 

0.134  X101 

256 

0.626  XlO"15 

0.119  xl0~14 

0.106  x  10~14 

0.426  X10~13 

0.262  XlO1 

512 

0.635  xl0“15 

0.149  XlO’14 

0.934  x  10-15 

0.426  X10~13 

0.520  xlO1 

Table  2:  Numerical  results  for  Example  1,  p  -  16. 


n 

£2(0) 

E°°(0) 

£2(0') 

E~W) 

t  (sec.) 

24 

48 

96 

0.764  XlO"10 
0.970  XlO'15 
0.851  X10"15 

in 

Brafrvffl 

Hi 

0.380  xlO0 
0.720  XlO0 
0.144  XlO1 

Table  3:  Numerical  results  for  Example  1,  p  =  24. 


Example  2.  The  purpose  of  this  example  is  to  demonstrate  the  performance  of  the  method 
when  the  coefficients  p,  q  of  the  equation  (6)  are  singular  at  the  ends  its  interval  of  definition, 
while  the  particular  solution  being  sought  is  smooth.  We  solve  the  Bessel  equation 

4>"(x)  +  ±-4>'{x)  +  ?—^=  0  (98) 

X  X i 

on  the  interval  [0.600]  with  the  boundary  conditions 

<t>{  0)  =  0, 

0(600)  =  1, 
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(99) 

(100) 


and  v  =  100.  The  difficulty  of  this  problem  is  due  to  the  fact  that  the  two  linearly  independent 
solutions  of  equation  (98)  are  Ju(x)  and  V^x)  (Bessel  functions  of  the  first  and  the  second 
kinds,  respectively).  As  is  well  known,  J„(x)  behaves  in  the  vicinity  of  zero  like  x",  while  V„(x) 
behaves  like  x-l/;  most  methods  have  trouble  finding  the  decaying  solution.  In  addition,  this 
is  a  fairly  large-scale  calculation,  since  the  interval  [0,600]  contains  almost  100  wavelengths  of 
the  solution  to  (98).  The  algorithm  has  been  applied  to  this  equation  with  p  =  16,20,  and  24. 
The  results  are  presented  in  Tables  4-6. 


T? 

E2(o) 

E°°(<p) 

£V) 

(<*>') 

t  (sec.) 

192 

0.945 

0.103 

0.101  xlO1 

0.172  xlO1 

0.194  xlO1 

384 

0.651 

0.901  xlO"1 

0.658 

0.108  xlO1 

0.394  XlO1 

768 

0.106  xlO-3 

0.151  xlO"4 

0.106  xlO'3 

0.177  XlO'3 

0.786  XlO1 

1536 

0.284  XlO-8 

0.406  xlO'9 

0.285  xlO'8 

0.478  xl0~8 

0.156  xlO2 

0.179  XlO'10 

0.265  xlO-11 

0.177  xlO'10 

0.310  xlO~10 

0.314  XlO2 

6144 

0.100  XlO"10 

0.138  x  10~n 

0.101  xlO~10 

0.167  xlO'10 

0.625  XlO2 

Table  4:  Numerical  results  for  Example  2,  p  =  16. 


n 

E2(o) 

Ex(0) 

E\<t>') 

£*(<?') 

240 

0.120  xlO1 

0.198 

0.111  xlO1 

0.232  xlO1 

0.300  xlO1 

480 

0.845  XlO"2 

0.118  xl0“2 

0.851  xlO-2 

0.140  XlO"1 

0.598  XlO1 

960 

0.684  XlO-7 

0.979  xl0~8 

0.687  xlO*7 

0.114  XlO"6 

0.119  XlO2 

1920 

0.205  xlO"11 

0.302  xlO-12 

0.202  xlO-11 

0.355  XlO"11 

0.239  XlO2 

3840 

0.229  XlO-10 

0.325  xlO"11 

0.231  xlO"10 

0.382  XlO"10 

0.480  XlO2 

Table  5:  Numerical  results  for  Example  2,  p  =  20. 


n 

E2(<t>) 

E°°(4>) 

E2m 

E™(<t>') 

t  (sec.) 

288 

0.889 

0.113 

0.942 

0.155  XlO1 

0.430  xlO1 

576 

0.765  XlO'4 

0.108  XlO'4 

0.770  xlO"4 

0.127  xlO'3 

0.866  XlO1 

1152 

0.206  xlO"10 

0.295  XlO'11 

0.207  xlO"10 

0.346  XlO"10 

0.176  XlO2 

2304 

0.356  XlO"11 

0.503  xlO'12 

0.356  xlO”11 

0.594  XlO"11 

0.343  xlO2 

4608 

0.627  xlO-11 

0.856  xlO-12 

0.644  xlO"11 

0.982  xlO"11 

0.688  XlO2 

Table  6:  Numerical  results  for  Example  2,  p  =  24. 


Remark  5.1.  Problems  like  the  preceding  one  are  frequently  encountered  in  the  modeling 
of  wave  phenomena  by  means  of  separation  of  variables,  and  were  the  original  motivation  for 
this  work. 
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n 

E\q) 

£°°(o) 

W) 

t  (sec.) 

240 

0.706 

0.913 

0.907 

0.675  XlO6 

0.242  xlO1 

256 

0.960  x  10-2 

0.468  xlO"1 

0.455  xlO'1 

0.338  XlO5 

0.256  xlO1 

272 

0.701  x  10-4 

0.445  XlO'3 

0.476  XlO-3 

0.265  xlO3 

0.278  XlO1 

28S 

0.835  xlO-7 

0.635  xl0~6 

0.761  XlO-6 

0.404 

0.288  xlO1 

304 

0.198  xlO"10 

0.141  XlO'9 

0.200  XlO"9 

0.147  Xl0~3 

0.308  xlO1 

320 

0.378  xlO"11 

0.233  xlO"10 

0.254  x  10~10 

0.289  XlO-4 

0.320  xlO1 

336 

0.610  xlO"11 

0.394  xlO~10 

0.321  XlO"10 

0.431  Xl0~4 

0.338  xlO1 

Table  7:  Numerical  results  for  Example  3,  p  =  16. 

Example  3.  We  solve  a  singular  perturbation  problem  of  the  form 

f  •*"(*)  -<t>'(x)  =  0,  (101) 

#-l)=l,  (102) 

0(1)  =  2,  (103) 

with  €  =  10~6.  The  solution  of  this  problem  has  an  extremely  sharp  boundary  layer  near 

the  right  end  of  the  interval  [  —  1,1],  causing  severe  numerical  difficulties  when  most  standard 

algorithms  are  used.  In  this  case,  we  construct  the  intervals  Bj71  =  [f>:,6,+i]  via  the  formula 

6i  =  —  1  , 

f>,  =  +  for  i  =  2,...,  M  ,  (104) 

&A/+1  =  1  , 

so  that  they  become  progressively  smaller  near  the  right  end  of  the  interval  [-1,1].  The  results 
of  this  experiment  are  presented  in  Table  7. 

Example  4.  Here,  we  solve  a  problem  of  the  type  which  arises  when  dealing  with  a 
frequency  domain  equation  for  the  vibrating  string. 

<t>"(x)  +  k2  •  <j>(x)  =  5  •  sin (k  ■  x),  (105) 

with  Diricmet  boundary  conditions 

4>(c)  =  sin (k  ■  c),  { 10G ) 

<p(d)  =  sin (k  ■  d)  .  (107) 

These  boundary  conditions  correspond  to  the  solution  <j>  =  sin(k  ■  x).  In  Tables  8-10,  we 
present  the  results  obtained  by  with  our  algorithm  for  c  =  -l,d  =  1,  and  k  =  630,  in  order  to 
demonstrate  the  performance  of  the  method  on  large-scale  oscillatory  problems  (the  string  is 
roughly  200  wavelengths  long). 
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n 

E\0) 

£°°(0) 

£'V) 

t  (sec.) 

800 

■ 

0.402  101 

0.353  101 

■ 

0.820  10' 

1600 

0.101  lO"3 

0.887  10~4 

0.165  102 

3200 

0.281  10"8 

0.242  10~8 

0.329  102 

6400 

EEESl 

0.187  lO"8 

0.158  10~8 

0.664  102 

Table  8:  Numerical  results  for  Example  4.  p  =  16. 


n 

£2(<p) 

£°°(0) 

EH#) 

E">W) 

t  (sec.) 

g 

0.218  10' 
0.839  10~4 
0.206  10~10 
0.411  lO"10 

0.379  10’ 
0.104  lO"3 
0.361  10-'° 
0.916  10-'° 

0.215  10' 
0.918  10"4 
0.181  10-'° 
0.353  lO"10 

0.241  104 
0.697  lO"1 
0.153  10"7 
0.352  lO"7 

0.894  10' 
0.179  102 
0.360  102 
0.721  102 

Table  9:  Numerical  results  for  Example  4,  p  =  24. 

Tl 

mtMsm 

E°°(<t> ) 

E2(<t>') 

E°°(<y ) 

t  (sec.) 

800 

1600 

3200 

6400 

0.223  10-' 
0.607  10~9 
0.194  10-'° 
0.118  10-9 

0.400  10"' 
0.837  lO"9 
0.252  10-'° 
0.172  lO*9 

0.235  10-' 
0.664  10"9 
0.168  10"'° 
0.112  lO"9 

EH 

0.164  102 
0.328  102 
0.675  102 
0.133  103 

Table  10:  Numerical  results  for  Example  4,  p  =  32. 


The  following  observations  can  be  made  from  Tables  1  -  10,  and  are  corroborated  by  our 
more  extensive  experiments. 

1.  The  practical  convergence  rate  of  the  method  is  consistent  with  the  theoretical  one.  For 
larger  p,  the  exact  numerical  verification  of  the  order  of  convergence  tends  to  be  difficult,  since 
the  precision  of  calculations  is  exhausted  before  the  behavior  of  the  scheme  becomes  asymptotic. 
However,  this  is  often  encountered  when  dealing  with  rapidly  convergent  algorithms. 

2.  For  small-scale  problems  (such  as  in  Example  1)  and  large  p,  the  algorithm  produces  essen¬ 
tially  exact  results  with  a  small  number  of  nodes.  For  large-scale  problems,  double  precision 
accuracy  is  achieved  at  approximately  20  nodes  per  wavelength  with  p  =  20,  at  12  nodes  per 
wavelength  with  p  =  24,  and  at  10  nodes  per  wavelength  with  p  =  32.  The  optimal  timings  are 
achieved  at  p  between  24  and  32  (provided  that  about  10  -  12  digits  of  accuracy  are  desired). 
There  seems  to  be  no  reason  for  using  the  scheme  with  p  <  16. 

3.  The  scheme  is  completely  indifferent  to  the  extreme  stiffness  near  the  left  end  of  the  interval 
[c,d]  of  equation  (98)  in  the  Example  2. 

4.  It  is  easy  to  use  the  algorithm  in  an  adaptive  manner,  as  demonstrated  in  Example  2.  where 
we  resolve  a  boundary  layer  of  relative  thickness  10-6,  without  encountering  any  numerical 
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difficulties.  However,  a  fully  adaptive  version  of  the  scheme  has  not  been  implemented.  The 
intervals  B™  in  Example  2  were  provided  by  the  calling  program  (as  opposed  to  having  been 
constructed  by  the  algorithm  itseif). 

5.  The  condition  number  of  a  Nystrom  discretization  of  a  second  kind  integral  equation  is 
asymptotically  bounded,  and  our  results  reflect  this  fact.  The  relatively  poor  accuracy  (10  - 
11  digits)  obtained  in  Examples  2,  3  and  4  is  due  to  the  ill-conditioning  of  the  original  ODE. 
as  opposed  to  that  of  the  numerical  scheme  used. 

VI.  Generalizations  and  Conclusions 

Generalizations:  The  results  of  this  paper  can  be  generalized  in  three  obvious  directions. 

1.  As  described  here,  the  algorithm  is  only  applicable  to  boundary  value  problems  for  a  single 
second  order  ODE.  It  can  be  generalized  to  systems  of  ODEs  of  arbitrary  dimension.  This 
generalization  is  quite  straightforward,  though  technically  somewhat  involved.  This  work  is 
currently  in  progress,  and  its  results  will  be  reported  at  a  later  date. 

2.  The  algorithm  of  this  paper  can  be  used  as  a  solver  for  the  linearized  problems  which 
arise  in  applying  Newton-type  methods  to  non-linear  boundary  value  problems  of  ordinary 
differential  equations.  This  would  involve  converting  the  differential  equation  being  solved  into 
a  non-linear  second  kind  integral  equation,  with  the  subsequent  application  of  an  iterative 
(Newton)  algorithm  to  the  latter.  In  the  context  of  non-linear  problems,  second  kind  integral 
equations  retain  their  usual  analytical  and  numerical  advantages  over  the  differential  equation 
formulation.  However,  our  numerical  experience  with  such  problems  is  extremely  limited. 

3.  Attempts  have  been  made  to  extend  the  results  of  this  paper  to  elliptic  partial  differential 
equations.  While  this  direction  of  research  appears  to  be  extremely  attractive  in  principle,  we 
have  not  been  able  to  produce  a  workable  algorithm  of  this  type. 

Conclusions:  An  algorithm  has  been  presented  for  the  solution  of  two-point  boundary 
value  problems  of  ordinary  differential  equations.  The  algorithm  is  based  on  reducing  the 
differential  equation  to  a  second  kind  integral  equation,  with  the  subsequent  solution  of  the 
latter  via  a  Nystrom  type  scheme.  It  has  CPU  time  requirements  proportional  to  N  p2,  where 
A  is  the  number  of  nodes  in  the  discretization  of  the  interval  of  definition  of  the  equation,  and 
p  is  the  desired  order  of  convergence  of  the  scheme.  The  method  does  not  involve  the  solution 
of  linear  systems  with  large  condition  numbers,  permits  the  use  of  schemes  with  extremely  high 
orders  of  convergence,  and  is  quite  insensitive  to  boundary  layers  or  to  end-point  singularities 
in  the  coefficients  of  the  differential  equation. 

Appendix  A. 

A  High  Order  Scheme  for  the  Solution  of  Equations  (82)  -  (84) 

In  this  appendix,  we  summarize  a  classical  approach  to  the  solution  of  integral  equations 
via  the  Nystrom  algorithm  based  on  Chebyshev  quadrature.  The  facts  used  in  this  appendix 
are  well-known,  and  can  be  found,  for  example,  in  [4,5,7]. 


23 


A.l.  Chebyshev  approximation. 

For  any  non-negative  integer  n,  the  Chebyshev  polynomial  Tn  of  degree 

n  is  defined  by  the 

formula 

Tn( cos  0)  =  cos(n$) . 

(108) 

Clearly,  |T„(x)|  <  1  for  x  £  [-1,1], 

T0{x)  =  1.  Ti(x)  =  x 

(109) 

and.  using  elementary  trigonometric  identities, 

Tn+i(x)  —  2xTn(x)  —  Tn—\(x)  for  n  >  1  . 

(110) 

It  is  easy  to  see  that,  for  n  >  1,  the  roots  I),,*2,  •  •  -,f£  of  Tn  are  real,  located  on  the  interval 
[-1, 1],  and  given  by  the  formula 


t'n  =  cos 


(2»-  l)ff 
2tj 


(111) 


The  Chebyshev  polynomials  constitute  an  orthonormal  basis  for  T2[-l,l]  with  respect  to  the 
inner  product 

U-g)r  =  J  ^f(t)  ■  g{t)  ■  {l  -  x2)~?  dt,  (112) 

Therefore,  any  function  /  £  C°(-l,  1]  can  be  represented  by  an  expansion 

/(*)  =  'T>(x)’  (113) 

i=0 

with  the  coefficients  a,  given  by 

a,  =  (/,T,)t-  (114) 

The  popularity  of  Chebyshev  expansions  as  a  numerical  tool  is  largely  a  consequence  of  the 
following  two  lemmas.  The  first  demonstrates  that  the  series  converges  rapidly  for  sufficiently 
smooth  functions,  while  the  second  shows  that  numerical  evaluation  of  the  coefficients  takes  a 
particularly  simple  form.  Proofs  may  be  found  in  [7], 

Lemma  A.l.  Suppose  that  n  and  k  are  natural  numbers,  and  that  f  £  Cfc[— 1,  l].  Suppose 
further  that  the  coefficients  ao,  aj,  •  •  ■  ,an  are  defined  by  formula  (114)-  Then  for  any  x  £  (—  1. 1]. 

|/(x)-]Ta.-r<(x)|  =  0(-7^T)  .  (115) 

n 

In  particular,  if  f  £  C°° ,  then  the  expansion  (113)  converges  to  f  superahebraically. 


V 
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A  well-known  property  of  the  coefficients  of  the  Chebyshev  series  is  that  they  can  be  ob¬ 
tained  from  the  cosine  transform  of  the  function  F(9)  =  f(cos9).  More  precisely,  a  simple 
change  of  variables  yields  the  result 

a0  =  (/. T0)t  =  ~  (  f(cos9)d9 
tr  Jo 

ak  =  {f-Tk)r  =  —  f  /(cos  9)  cos k9d0  for  k  >  0.  (110) 

JT  Jo 

Lemma  A. 2.  Suppose  that  n  and  k  arc  natural  numbers,  that  f  E  C*[-l,l],  and  that 
the  coefficients  a,  are  defined  by  (114)-  Suppose  further  that  the  vector  fn  =  (/„,/„.•••,/„) 
consists  of  the  function  values  at  the  roots  ofTn(x),  namely 

/;  =  /K),  i  =  1,2....,  A.  (117) 

Let  q  =  ( Oo , Q i , •  •  ■ , Qn_i)  be  given  by 

a=Cn(fn),  (11 S) 

where  Cn  denotes  the  discrete  cosine  transform  of  dimension  n.  Then 

\a,  -  a,\  =  0(±)  • 
nK 


A. 2.  Chebyshev  quadrature. 

Lemmas  A.l  and  A.2  provide  a  tool  for  the  construction  of  highly  accurate  interpolation 
schemes.  The  following  two  lemmas  use  Chebyshev  expansions  as  an  apparatus  for  the  numer¬ 
ical  evaluation  of  indefinite  integrals. 

Lemma  A. 3.  Suppose  that  f  :  [-1,1]  — *  R  is  given  by  the  finite  Chebyt  ci  ..,les 

n- 1 

/(*)  =  -r,(x)  .  (119) 

i=0 

Then  the  indefinite  integral  of  f  has  a  series  expansion  of  the  form 

Fi(x)=  [T  f(t)  df  =  f>  r,(x)  .  (120) 

The  coefficients  are  given  by 


P>  =  YTi  *  (Ci-1  a.-i  “Ci+1  al+1)  for  t  >  1  , 

00  =  2-£(-l  r1  -/J,.  (121) 

1=1 
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where  cq  =  2,  c,  =  1  for  i  >  0,  and  a,  is  assumed  to  be  zero  outside  the  range  i  =  0, 1 , . . . ,  n  —  1 . 


Definition  A.l  Let  a  —  (qo,  ct\.  ■  ■  ■ ,  Q„-i)  and  let  0  =  (0o,  3 1,  Then  the  linear 

mapping  Sn  :  Rn  — •  Rn+1  defined  by  the  formula 

Sn(a)  =  3 


is  referred  to  as  the  spectral  integration  operator. 

Lemma  A.4.  (Clenshaw  and  Curtis  [2])  Suppose  that  f  6  [ — 1 , 1]  for  some  k  >  1,  and 

that  the  vector  fn  =  (/,},/,?,. .  -  G  Rn  consists  of  the  function  values  at  the  roots  ofTn(x). 

/;  =  /kk  <  =  a'. 

Suppose  further  that  Fi  :  [-1, 1]  — *  R  is  the  indefinite  integral 

Ft(x)  =  j\mdt 

and  that  the  vector  Fn  =  (F,},  F%, . . F")  £  Rn  is  defined  by  the  formula 

F^Ffit'J.  (122) 


Then 

||F"  -  C-1  o  S’1  oCn  (Dlloc  =  0(^y )  •  (123) 

Furthermore,  all  elements  of  the  matrix  C'1  o  Sn  o  Cn  are  strictly  positive. 

The  following  theorem  provides  a  tri  via]  extension  of  the  above  lemma  to  the  case  of  intervals 
of  length  other  than  2. 

Theorem  A.l.  Suppose  that  [a,i]  C  R  is  an  interval  of  non-zero  length,  n,k  is  a  pair  of 
natural  numbers,  f  £  C*[a,  6],  and  the  finite  sequence  t*,  t^,  •  •  •  ,  r”  is  defined  by  the  formula 


cos 


f(2T  -  1)7T 


2n 


+ 


fb  +  a 

V  2 


Suppose  further  that  the  function  Fi  :  [a,  6]  — *  R  is  defined  by  the  formula 

f,(x)=  r  f(t)dt. 

Ja 


(124) 


(12'. ) 


Finally,  suppose  that  the  vectors  fn  -  {f\,  f%, . . . ,  f" )  and  Fn  =  (F,J,  F*, . . F" )  are  defined 

h 


Then 


(128) 


II  Fn~ 


b  —  a 


Cn  0  Sn  a  Cn  (/")|!  =  0(- 


U 


In  the  solution  of  the  integral  equation  (8),  we  will  also  require  the  evaluation  of  indefinite 
integrals  of  the  form 

Fr(*)  =  £f(t)dt. 

For  this  purpose,  we  have  the  following  theorem. 

Theorem  A.2.  Suppose  that  [a,  6]  C  R  is  an  interval  of  non-zero  length,  n,k  is  a  pair  of 
natural  numbers,  f  E  Ck[a,b ],  and  the  finite  sequence  rf,  r%,  •  ■  • ,  r”  is  defined  by  the  formula 


(2 i  -  1)t‘ 
2n  . 


Suppose  further  that  the  function  Fr  :  [a.  6]  — >  R  is  defined  by  the  formula 

Fr(x)  =  [bf(t)dt. 

J  X 


(129) 


(130) 


Finally,  suppose  that  the  vectors  fn  =  (/£,/*,...,/£)  and  Fn  =  (iqj,  F%, . . . ,  F„  )  are  defined 

by 

fk  =  f«)  (131) 

and 

K  =  (132) 

Then 

\\Fn  -  ~  ■  C-1  oSnoCn  (Dll  =  0(- 1T) ,  (133) 

where 

Sn  =  T  o  Sn  o  T  (134) 

and  T  :  R"  — *  Rn  is  the  transposition  operator  defined  by 


T’l.n-i+l  ~  li 

Tn-i+\,i  =  1, 


T 


0 


otherwise. 


(135) 


Definition  A. 2.  The  matrix 


on  _ 
^a,b  ~ 


•  C~l  o  Sn  o  C„ 


(136) 
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is  known  as  the  left  spectral  integration  matrix  of  order  n  on  the  interval  [a,  6].  The  matrix 

Sib  =  b~~  ■  C~X  oSnoCn  (137) 

is  known  as  the  right  spectral  integration  matrix  of  order  n  on  the  interval  [a,  6]. 

According  to  Theorems  A.l  and  A. 2,  S£t  converts  the  values  of  a  function  /  :  [a,&]  —  R  at 
the  Chebyshev  nodes  on  the  interval  [a,  6]  into  (approximate)  values  of  the  indefinite  integral 

r  m* 

J  a 

at  these  same  nodes.  Similarly,  (,  converts  the  values  of  /  into  (approximate)  values  of  the 
indefinite  integral 

l  M  dt  • 

Furthermore,  the  order  of  convergence  of  the  discrete  approximations  is  equal  to  the  number  of 
continous  derivatives  /  possesses.  Finally,  S^b  and  can  be  applied  to  an  arbitrary  vector 
for  a  cost  proportional  to  n  •  log(n),  since  it  is  a  product  of  a  diagonal  matrix  ( Sn  or  S")  and 
two  cosine  transforms.  The  latter  fact,  however,  is  irrelevant  for  the  purposes  of  this  paper. 

A.3.  The  Nystrom  algorithm  for  the  solution  of  second  kind  integral  equations. 

The  Nystrom  algorithm  associated  with  a  family  of  n-point  quadrature  formulae  = 
{zk,w,}  replaces  the  integral  equation 

4>(z)+  f  K(z,t)  ■  4>{t)dt  =  g(z)  (138) 

Jo 

with  the  system  of  linear  algebraic  equations 

n 

4>,  +  Ylw'r  K(zt,Zj)  ■  <t>}  =  g(zx)  t  =  1, . . ., n.  (139) 

i= i 

We  denote  the  matrix  of  the  discretized  system  (139)  by  Bn ,  and  view  the  solution  <p\ , . . . ,  4>n 
of  (139)  as  an  approximation  to  the  solution  <f>  of  (138)  at  the  nodes  zi,...,zn.  If  (139) 
has  a  unique  solution,  then  for  a  wide  class  of  quadrature  formulae  rjn  and  sufficiently  large 
n,  the  system  (138)  also  has  a  unique  solution.  Furthermore,  under  fairly  broad  assumptions, 
the  convergence  rate  of  the  Nystrom  algorithm  is  the  same  as  the  convergence  rate  of  the 
quadrature  formula  it  is  based  on  (see,  for  example,  [1]). 

A.4.  A  high  order  Nystrom  scheme  for  the  solution  of  equation  (8). 

In  this  section,  we  describe  a  particular  version  of  the  Nystrom  scheme  for  the  solution 
of  equations  of  the  form  (8).  The  scheme  to  be  described  is  based  on  Chebyshev  quadrature 
and  has  a  convergence  rate  of  the  order  k  —  1,  where  k  is  the  number  of  continous  derivatives 
possessed  by  the  solution  o  of  equation  (8). 
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Given  an  interval  [<7,c]  and  the  Chebyshev  nodes  tj,T2,  •  • -,rn  on  it,  we  will  define  the 
diagonal  operators  P,Q,Ul,UT ,  Ul\UT',Vl,VT  :  Rn  Rn  via  the  formulae 


Pi.i 

= 

P(a), 

Qi.i 

= 

9(a). 

= 

Ui(r,). 

uh 

= 

Ur(T>), 

vl 

= 

«j(a), 

K 

= 

<(a). 

VU 

Vi(t,), 

vr 

»,* 

= 

Vr  (A). 

Finally,  we  will  define  the  operator  An  : 

Rn 

— »  Rn  by  the  formula 

An  =  I  +  Po^o^oFU^'o^or), 

+  +  (141) 

and  the  vector  fn  £  Rn  by  the  formula 


f?=f(r,).  (142) 

Observation  A.l.  The  mapping  A"  defined  by  the  expression  (141)  can  be  viewed  as  an 
approximation  to  the  operator  defined  by  the  expression  (81).  This  is  obvious  from  (11). 
(12),  Theorem  A.l  and  Theorem  A.2. 

The  following  theorem  provides  an  exact  statement  of  the  above  observation.  It  is  used 
in  Section  IV  of  this  paper  to  construct  highly  accurate  approximations  to  the  solutions  of 
equations  (82),  (83)  and  (84)  on  elementary  subintervals  B™ ,  and  is  the  principal  purpose  of 
this  Appendix. 

Theorem  A. 2.  Suppose  that  the  equation  (7)  is  chosen  in  such  a  manner  that  tv.iy  £ 
C*[a,c]  for  some  k  >  1.  Suppose  further  that  the  equation  (82)  (or  (83)  or  (84))  has  a  unique 
solution  a  such  that  a  €  C*[a.c].  Then  for  all  sufficiently  large  n,  the  equation 

An(<rn)  =  fn  (143) 


has  a  unique  solution  an,  and 


Ik?  -  *(%)lloo  =  o(-xnO  • 

n K  1 


(144) 
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